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Abstract. We describe three algorithms for computer-aided symbolic multi-loop calculations 
that facilitated some recent novel results. First, we discuss an algorithm to derive the canonical 
form of an arbitrary Feynman integral in order to facilitate their identification. Second, we 
present a practical solution to the problem of multi-loop analytical tensor reduction. Finally, 
we discuss the partial fractioning of polynomials with external linear relations between the 
variables. All algorithms have been tested and used in real calculations. 



1. Introduction 

A higher-order calculation is a multi-stage process, with details strongly dependent on the 
physics problem. There is no generic one-size-fits-all approach, but usually one first generates 
the diagrams, then performs Dirac algebra and/or projections on scalar integrals, and then 
manipulates large expressions that depend on scalar products of loop and external momenta 
until they can be reduced to known integrals and computed. Each calculation has its own 
"key" stage which requires the most efforts. In many cases, those are the algebraic reduction of 
integrals (with integration-by-parts identities [H [21 [3] ) to a small subset of "master" integrals, 
and the evaluation of the latter. Those problems have been the focus of active research for 
many years and gave rise to a plethora of various methods and algorithms [4J. However, as the 
three-, four- and higher- loop problems become the norm in the LHC era, even the traditionally 
"simple" steps become impossible to perform manually or with the ad-hoc methods. 

In this contribution we present the three algorithms that automate the tasks usually 
considered "routine" and not worth mentioning. First, we discuss the classification of integrals 
as instances of different "topologies", or families of integrals that differ by the exponents of 
the (fixed) denominator factors. When the number of such families reaches a few hundred, 
automation becomes a necessity. 

Second, we consider the tensor reduction, that may e.g. help re-arrange numerators of 
factorizable topologies in order to perform integration independently. In certain cases, involving 
complicated asymptotic expansions, this stage may dominate the total computation time of a 
diagram. In addition, we provide a practical solution for the situations where general formulas 
are not known. Finally, we present a (relatively straightforward) re-formulation of the partial 
fractioning problem and solve it using Grobner bases. 

The primary use of these methods is to generate the code for the computer algebra system 
that processes the actual diagrams (e.g. FORM ^), and thus the runtime efficiency is not the 



highest priority. Nevertheless, we find that our solutions are rather efficient and do not require 
enormous times for real-life problems. 



2. Identification of Feynman integrals 

In multi-loop calculations, one often has to solve a problem of identifying individual Feynman 
integrals or deciding whether an integral belongs to a given family ("topology") which helps 
reduce the number of integrals to compute. In general, in order to say whether two integrals are 
equal one has to compute both. However, for some integrals it is possible to elaborate a simple 
transformation of loop momenta and establish the identity of integrands. In simple cases, such 
transformations can be derived manually just by looking at graphs, but in multi-loop multi-scale 
problems with non-trivial kinematic constraints this is a daunting task. 
Let us consider two forward-scattering integrals (Fig. [1]) : 
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with the common additional constraints pf = P2 = 0, {pi +^2)^ = — 1- 

The integrands of Ii and I2 obviously cannot be related by only renaming symbols. However, 
the value of e.g. Ii is invariant with respect to a shift of loop momenta, e.g. ki k[ = ki+pi, 
or a permutation of denominator factors, e.g. Di o D2, etc. There may exist a transformation 
Di{K) = Ep(^i^{Q) with K = MQ and p{i) is some permutation of indices 1, 2, 5, such that 
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A and B are the 2x2 matrices, \A\ = 10 , and / is a 2 x 2 identity matrix. 

One practical way to compare two topologies is to compare graphs. In Fig(T]we present the 
graphs corresponding to Ii and I2, and one can see that the graphs do not match. However, 
some integrals may have no corresponding graphs (for example, integrals originating from some 
effective theories may have linear denominators, e.g. 2piki + iO, which cannot be interpreted as 
graph lines), and some integrals may have multiple corresponding graphs (as in the given case). 
In addition, (sub)graph isomorphism is itself a computationally hard problem^. 

^ or in general \A\ = d ^ 0, and integrals then are equal modulo corresponding Jacobian. 

^ Subgraph isomorphism problem is NP-complete, while graph isomorphism belongs to the NP complexity class. 




Figure 1. Graphs corresponding to integrals Ii and I2. 



Instead, in order to compare integrals, one could use alpha-representation [6j, which is 
explicitly covariant, very closely related to the definition of dimensionally regularized Feynman 
integrals and can be derived for any set of quadratic denominators. The structure of any 
Feynman integral is encoded by the two homogeneous polynomials, U and F which do not depend 
on the exponents of denominator factors and the dimensionality of space D. In particular, we 
have: 



with constants C, a and b (irrelevant here) depending only on D. 

The only freedom that is left when comparing {Ui,Fi} with {U2, F2} is a permutation of lines. 



mundane properties of polynomials common to all permutations, such as the total number of 
terms or the number of terms proportional to m^. 

An obvious way to find p{i) when comparing two integrals would be to try all 5! = 120 
permutations of Xj. Since usually one has more than two candidates to compare, it is beneficial to 
introduce a "canonical" ordering, maximizing some metric over all permutations. The canonical 
ordering of lines should then be derived only once per integral, and the identity of canonical 
alpha-representations would provide the definitive answer. 

In practice, we found that it is sufficient to build a metric in the space of not pairs {U, F} 
but products UF. One example of a suitable metric is given by the following rules (we assume 
that a unique ordering of coefficients is available, as it is in any computer algebra system): 

(i) Turn the polynomial of n variables with m terms into a matrix with rows corresponding to 
monomials. The first column contains coefficients, and the subsequent columns contain the 
(non -negative integer) exponents of variables xi, ...,Xj7,. 

(ii) Make n copies of the table. In the i-th copy, exchange the second column (corresponding 
to xi) with the i-th column (originally corresponding to Xj). 

(iii) In all copies, sort rows lexicographically by the first two columns (i.e. compare only the 
first two entries in each row). 

(iv) Extract the second column from each copy (as vectors of length m), and determine the 
lexicographically largest vector, comparing all m elements. 

(v) In the table copies with the maximized second column, continue recursively: produce n — 1 
copies-of-copies, in each select a different third column, sort by the three first entries, find 
the maximum third column, discard non-maximal entries, etc. 

(vi) The permutations of columns in the copies maximizing all columns (there can be a few due 
to symmetries) can be taken as the "canonical permutations" of Xj. 

While there is a theoretical possibility of combinatorial growth in this strategy, we 
have found this approach quite fast and practical to at least five-loop integrals. In our 
example, after the canonical re-ordering, one easily establishes that f7iFi(x5, X2, X3, xi, X4) = 
^^2-^2(2:5, xi, X4, X2, X3), i.e. Ii and I2 are indeed identical. After this identification it is easy to 
find the relation between the loop momenta. The graphs, however, may differ as shown in Figj2j 
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Ui = U2 = X5(X3 + X4) + (Xi + X2)(X3 X4 X5), 

Fi = m^x^Ui - X1X4X5, F2 = m^x^Ui - X2X3X5. 




3. Analytic tensor reduction 

A traditional approach to computing multi-loop integrals with tensor structures in numerators 
is to decompose the result into all tensors allowed by the symmetries and then determine the 



Figure 2. Graphs corresponding to integrals Ii and I2 after identical loop re-parametrization. 



coefficients of the decomposition. Let us consider a simple one-loop example: 

/ A:'^^A;^2^'^3^M4 . = {k''\..k>'^) = Y CiTl'^-^\ (5) 

The tensors xj^^'"^'^ cannot depend on any external momenta, thus they may only be composed 
of metric tensors. There are only three unique products: Ti = g^^^^^'^g^^s,^^'l ^ j^^ — g^^l^^3g^^2^^4, ^ g^j^^ 

The coefficients Q can be found by multiplication of the original equation Eq. [7]with Ti, T2, 
and T3 and solving the system of equations: 

Ti^,,...^,{k^' ...k^') = {{k'f) = CiD^ + C2D + CsD, (6) 
T2^,...>.Ak^'-k^') = {{k"?) = CiD + C2D^ + C3D, 
T3^,...>.Ak^'-k^') = {{k^) = CiD + C2D + CsD^ 

The equality of all products on the left hand side implies that Ci = C2 = C3 = C, thus 
reducing the number of independent equations to one, and we 

findC7= [Z)(D + 2)]-^((fe2)2). 

In this case, it is not difficult to generalize the formula to a general number of indices: the 
result will be represented in terms of symmeterized product of metric tensors, and the coefficients 
will contain a number of gamma- functions. In a higher number of loops, however, there exist 
only a few general formulas. The two- loop generalization of EqjTl first given in [7] involves 
Gegenbauer polynomials, traceless tensors, and is much more computationally expensive due to 
multiple nested summations. Tensor reduction may also be derived in terms of dimension shifts 
and differential operators [8J, but that is not very convenient in real calculations. 

Instead, we suggest to directly use the method as outlined above for the fixed number of 
indices. For illustration, let us consider a two-loop propagator-type integral with loop momenta 
ki and k2, and the external momentum p, where the formulas of Davydychev-Tausk still apply 
(although with some effort), and limit ourselves with six open indices in the numerator. 

In the left-hand side we can have seven distinct distributions of indices over the loop momenta 
(e.g., {ki^ki^ki^ki'^k2^k2'^)). On the right-hand side we may build 76 possible tensors out of 
metric tensors and the components of the external momentum p (e.g., gl^^^^^ p^^3 gfi4fi5 pt^e^ _ 

For each of the seven LHS combinations, the coefficients in front of the RHS tensors have 
to be determined separately. Naively, in every case one would need to solve 76 equations with 
76 variables. However, the example above gives us a hint: one may reduce the number of 
independent variables by exploiting the symmetries. If the products of an LHS tensor with the 
two RHS tensors coincide, so will the corresponding coefficients. In this case, one has at most 
10 independent coefficients to determine (each time from 76 equations). 

This procedure does not depend on the actual denominators of the integral: we only derive 
a decomposition of the numerators. Thus, it is possible to pre-compute the tables with the 
reductions for a sufficiently large number of indices and re-use them in different computations. 



Figure 3. Original (linearly dependent, (a)) and fraction-decomposed ((b) and (c)) topologies. 



One only needs a practical way to solve the systems of linear equations, where coefficients depend 
on dimensionality D and possibly kinematic invariants. In our setup, we reuse the components 
from the Laporta algorithm that perform Gauss reduction. 

We have checked that this approach is indeed very practical: for example, our implementation 
builds a table for all 4-loop propagator-type integrals with up to 6 indices in the numerator in 
less than one hour with a regular PC. 

4. Partial fractioning 

Normally after the reduction to scalar integrals one has a long polynomial with terms that differ 
in exponents of denominator factors, e.g. 

diagram = DjDi + SSD^^D^^ - 45/16L>iL»^^ + (e + 2)D2 + (7) 
Di = /c^ + m^, D2 = k'^, i.e. linear dependence: Di — D2 — m? = 0. 

(the corresponding graph is in Fig. [3] (a)). In what follows we will discuss general transformations 
of monomials DfD2, where a and b are positive or negative integer numbers. 

The condition that Di and D2 are dependent allows us to simplify such monomials. For 
example, we may use the "classical" fraction decomposition relation 



D1D2 m'^D2 m?Di ^ ' 

repeatedly to reduce any monomial with negative exponents a and 5 to a linear combination of 
terms where at least one of a or 6 is non-negative. 

The form where one of D\ and D2 is absent is preferable for the further computations. In 
this case, the goal of the fraction decomposition is to identically re-arrange Eq. [7] so that each 
term in the final expression could be identified with one of the "simpler" topologies in Fig. EJ^b) 
and (c). 

However, the rule Eq. [8] alone does not achieve this goal. In particular, it does not apply 
when 6 > 0, a < and one has to apply the direct decomposition rule D2 D\ — rv? . In the 
symmetric case 6 < 0, a > it could be convenient to use D\ D2 + m^, but applying both 
those rules to monomials with a > and 6 > will not terminate. After the inspection of all 
relevant cases one can find that that only a system of three rules 

{DiD2r^ ^ {m^y^D^^ - D^^') , D2^Di-m^, Di/ D2 ^ / D2 + I (9) 

provides the complete decomposition. 

The example above is relatively straightforward and can be easily done manually. However, 
in practice one may have a much more complicated case for linearly dependent factors. In Fig. U] 
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Figure 4. Linearly dependent topology with two relations between denominator factors. 

is a two-loop example of a topology in special kinematics, where = — and each massive 
line has mass m. Here one has two relations between the seven factors: 

2Di -D^-D5- = 0, (10) 
2D2 - 2D3 + D4-D5 + 2L>6 - '^Dr + ^n? = 0. 

Let us formulate the problem in a slightly different way. First, let us get rid of the negative 
powers of Di by introducing variables Yi = D~^. Second, for all monomials D^^ ...D2^Y^^ ...Yj^'^ , 
dj > 0, let us introduce an ordering that incorporates our understanding of what is "simpler". 
In particular, it is beneficial to use linear weighting of exponent vectors (ai, ...,014). Given two 
such vectors a and b, we decide which of them is the "largest" by lexicographically comparing 
the products Ma and Mb, where M is some square 14 x 14 matrix. One reasonable choice is 
f 1 if i > j 

~ { n ■ "^^^ sum of all exponents is here the primary criterion, which agrees 



0, otherwise 

with the intuitive definition that the monomials with fewer non-zero exponent are "simpler" . 

Given this ordering, the problem is to re-arrange polynomial conditions Eq. [10] and the 
additional relations DiYi — 1 = into a set of equivalent relations that would unambiguously 
reduce any given monomial to the "simplest" form (and would not lead to loops during the 
repeated application). 

As formulated above, exactly this problem is solved by the so-called Buchberger algorithm, 
that generates a set of polynomials known as the Grobner basis. We then need to interpret each 
element of this basis (a polynomial p = 0) as a rule to substitute its most "complex" monomial 
(according to ordering M) with the remaining terms. Quite conveniently, Mathematica has 
a function GroebnerBasis [ . . .] that is rather efficient and has an option MonomialOrder to 
select the proper weight matrix M. In this particular case, this function produces a Grobner 
basis consisting of 14 polynomials that would be very difficult to derive manually. This output 
can be directly translated to FORM code and used to decompose expressions of any complexity. 
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